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Abstract 

Recent years have seen an increased level of interest in pricing equity options under a stochastic 
volatility model such as the Heston model. Often, simulating a Heston model is difficult, as a standard 
finite difference scheme may lead to significant bias in the simulation result. Reducing the bias to an 
acceptable level is not only challenging but computationally demanding. In this paper we address 
this issue by providing an alternative simulation strategy - one that systematically decreases the 
bias in the simulation. Additionally, our methodology is adaptive and achieves the reduction in bias 
with "near" minimum computational effort. We illustrate this feature with a numerical example. 
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1 Introduction 

Under the standard Black-Scholcs framework, the asset price dynamics is given by the lognormal 
model 

J Q 

—±=fidt + adW t (1.1) 

where the drift parameter \i and the volatility a are considered constant (or at best piecewise con- 
stant). The popularity of this model lies in its convenience and simplicity; however, these features 
come at a price. The standard Black-Scholcs model is not able to capture the volatility smile observed 
in the trading market. The Heston model assumes that the variance V = a 2 is itself a stochastic 
process - more specifically a square-root diffusion model of the CIR (Cox-Ingersoll-Ross) type (see 
[TTj . [5]). It further allows the variance, V, to be correlated with the stock price S, thereby captur- 
ing the volatility smile. Furthermore, the Heston model provides a closed-form solution for pricing 
European Options, allowing one to fit the model to observed option prices. The Heston model is 
given by the coupled SDE (stochastic differential equations): 



(idt + y/VtdWf (1.2) 



dfn 
S t 

d\'. = k{9 -V t )dt + a V y/v t dWy (1.3) 



where the variance V is modelled by a square-root diffusion process with parameters k which is 
the speed at which the process mean reverts to the long term variance 9, and oy, the volatility 
of the variance. We denote by p the instantaneous correlation between the two noise processes: 
d{W s ,W v ) =pdt. 
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By Cholesky factorization, one can rewrite equations (|1.2|l - (|1.3|l as: 



dSt 
S t 
dV t 



= pdt 
= k(6- 



VVt 
V t )dt 



pdW t {1} + y/l - fPdWt 



(2) 



with independent Brownian motions dW^ and dW\ L) . There is a wide range of literature on the sim- 
ulation of this model. Some of these are based on Euler discretisation; others on the improved finite 
difference approximations such as higher-order Milstein schemes and Predictor-Corrector methods; 
and still others based on distributionally exact ([5]) or approximate ([3]) simulation methods. In our 
study we focus mainly on the Exact simulation method as proposed by Broadie and Kaya [5] but we 
will avoid the numerical inversion of the Laplace Transform which seems to be the most time con- 
suming part of their simulation technique. Doing so, will mean introducing bias in the simulation. 
We will address this problem by providing an adaptive strategy that systematically controls the 
bias in the simulation. Furthermore, this is achieved with minimal computational overhead. In most 
cases, our method should prove faster than that of 0, as there is no costlier inversion to perform. 

Here is the layout of the rest of the paper. In Section [2j we recall the general algorithm (cf. Broadie 
and Kaya [5]) for simulating the Heston model, to emphasise the basic difficulty: the simulation of 
the integral of the variance process. In Section [3] we focus on the dynamics of the variance process 
and study its properties by transforming it into a canonical process. We also describe methods for 
simulating the latter process. In Section[3]we describe a method for simulating bridges corresponding 
to the latter process, and their use in simulating the integral of the variance process. In Section [5] 
we relate the integral of the variance process to a weighted integral of the canonical process, and 
we also calculate some moments of the latter for the canonical bridge process. In Sections [BJ 
we explore adaptive computation of the integral. In Section [5J we describe the results of some 
numerical experiments on the accuracy and efficiency of our adaptive algorithm. Finally, we provide 
concluding remarks and future research ideas in Section [9] Lengthy or highly technical proofs are 
given in Appendices A-C. 



<T V \/V t dWi 

(2) 



(i) 



(1.4) 
(1.5) 



2 Simulating the Stock Price 

All the details of the simulation algorithm for the equity price governed by the stochastic differential 
equation in (|1.4[> . are given in \5\. We reproduce them here for the sake of completeness and to 
emphasise the role played by the integral of the variance process, V. Integrating (| 1 .4[) between two 
dates of interest (e.g., coupon/reset dates), tj-\ and tj, we obtain 



Sj = Sj-i exp 



tj Ctrj 

V s ds + p 



^dW^ 



where Atj — tj — tj-\. Similarly, integrating (ll.5[) we obtain 



ftj rtn 

Vj = Vj-i + nOAtj — K I V s ds + oy 



fVsdW? 



(2.1) 



In the next section we will see that it is not complicated to simulate Vj in a manner that is not 
based on (|2.ip . We can easily use (|2.1[) to obtain: 



^VsdW^ = — [ AVj - nOAtj + n 



V s ds 



(2.2) 
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where AVj — V 3 ■ — Vj~\- The only component left is ^^/VsdWs 2 ^; but since V s is independent of 



(2) 

the Brownian increments dWs by construction, then given V s , < s < t 



j ■ 



ifv s dWi 2) ~ N 0, / F s ds 



Using this result we have: 



Sj = Sj^i exp 



1 



liAtj - - / V a ds - 




(l-p3) 



V, (is 



(2.3) 



where, given 14, < s < tj, Zj ~ AT(0, 1) and is conditionally independent of W^ 1 ). Based on (|2.2[) 
and (I2.3[) . it is clear that simulation of the Heston model is straightforward except for the simulation 
of the integral f* 3 ^V s ds. In |5 , it is done in an unbiased manner by using the distribution of this 
integral. The distribution is not available in closed form and is computationally intensive to compute. 
Our approach is to use a random quadrature with an adaptive control of bias. 



3 Squared Bessel Process and its Simulation 



Our first task is to be able to simulate equation (jl.5p exactly. In order to do so, we cite the following 
result from |12) : 

Theorem 3.1 Consider the one-dimensional diffusion process with laws {^Q d x ,x > 0} (x is the 
initial point; d > and f3 are fixed parameters, the former called the dimension), with infinitesimal 
generator 

2xD 2 + (2fix + d)D. 
Then for all real (3^0, is the Q^-law of the process, 



exp(2/3t)X 



1 - exp(-2/3i)\ 
W ) 



where Q d is the distribution of the d-dimensional squared Bessel Process denoted by BESQ 1 with 
inhnitesimal generator, 



2xD 2 + dD. 



Corollary 3.2 The space-time transformation 

V t = exp(-Kt)X T (f), V Q = x , (3.1) 
r(t) = -r[eM^)-l], T O =t Q =0 (3.2) 
transforms the square-root diffusion process in il.5\) to the X-dimensional squared Bessel process: 

dX u = Xdu + 2y r X\ l dW u , X{0) = x = V (3.3) 

where A = ^, 
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Proof A direct proof of this result is deferred to Appendix A. □ 

The squared Bessel process in p. 31) will also be referred to as a squared Bessel process of order 
v = A/2 — 1. From [4] we know that the boundary point is strongly reflective when — 1 < v < 
and is "entrance-not-exit" when v > 0. Following [TB] and [5] we simulate the squared Bessel process 
over any interval (tj,t,--|_i] with x(tj) known (already simulated for j > 0), using the randomized 
Gamma distribution of the first kind, G(v + rj + 1, 2At), where At = t j+ i — Tj and r\ is sampled 
from the Poisson distribution P(fj,),/j, = x(tj)/2At: 

X(t j+ i) ~ G(u + rj + 1, 2Ar) with r? ~ P 

Once we obtain the simulated values of X on the set of points Tj = T(tj),j = 1, 2, . . . , N, we obtain 
the corresponding Vj values using the transformation (|3.ip . 




4 Simulation of Bessel Bridge Process 

As already emphasised, the most important and difficult piece of our algorithm is the simulation of 
the integral J t J+1 V(s) ds. We propose to do this by recursively applying a Bessel bridge simulation, 
to fill in the intermediate points in the interval, [tj,tj+i], for a random quadrature. Methods for 
selecting the intermediate points are described later in the paper. In this section we describe the 
method of generation of a Bessel bridge process, and its application to the quadrature. 

Suppose that on any generic interval [tl,tr] we need to insert another point tm'- tl < Tm < tr. 
The corresponding BESQ X value xm is simulated using the randomized Gamma distribution of the 
second kind, G(-,-) (see e.g., [16], [8]): 

At \ 



X(r M )~G [v + m + 2r? 2 + l, 



2At l At r J 



where At = tr — tl, Atl = tm — tl , Atr = tr — tm , and where rji is sampled from a Poisson 
distribution, P(-), and 772 is sampled from a Bessel distribution, £>(•, •): 



Tjl ~ P 



1 



2At 



Atr At l 
At~ L XL + At~r XR 



„, y/X L XR 



Remark 4.1 For simulating the squared Bessel process and the Bessel bridge process we need 
efficient Poisson and Bessel random variate generators. For generating Poisson variates we refer to 
0, and [BJ- For generating Bessel variates we refer to Section 2 in and especially to the 
comments appearing in Section 4 in jl&j . 

Once we have generated intermediate values on K — 1 intermediate point^l ti,T2, tr_i [tq cor- 
responds to the left endpoint and tr corresponds to the right endpoint of the integration interval, 



1. These points may be chosen directly in the r-space, or in the t-space and mapped to r-space using the one-one 
mapping, A3. 21 . 
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[tj, tj+i]), the integral of the Variance process V may be approximated by its (conditional) expec- 
tations over the subintervals. Explicitly, with Sk = i" _1 (rfc), 



V(s) ds : 



K 

E 

fe=i 



E 



V(s) ds 



X(Tk-i) = xfa-i) ,X(T k ) = x(r k ) 



(4.1) 



In the next section, we develop a closed form formula for the conditional expectations on the right- 
hand side of (14. ip . For the sake of brevity, when the time interval is fixed, we sometimes denote the 
conditional expectation operator in (|4.1[) as E I[Ifl , where L and R signify the "frozen" left- and 
right-hand endpoints of the interval, respectively. 

We now have the strategy for simulating f t i+1 V(s) ds, for some given endpoints tj and except 
for the choice of K and the intermediate points. These topics are discussed in Section [5] 



5 Theoretical Results for Integrals of V, X 

In this section, we collect some theoretical results on the integrals of V and X. The proofs of some 
of the results are deferred to the Appendices. Our first result expresses the variance integral in the 
(V, £) space in terms of the one in (X, r) space. We work on a general interval, [£/, , tp>] which, in 
application, could be [tj,£j+i] or any of its subintervals upon refinement of a partition of [£j,£j+i]. 
Set 

tl = r(t L ), t r = t(£r). 



Lemma 5.1 On any interval [£l,£.r], the integral of the variance V is described in terms of the 
canonical BESQ X process X as: 

V t dt = / X(u)w a (u)du (5.1) 

wq(u) := <2q[1 + cw]~ ; ao := ,c:=4K(7y . (5-2) 



Proof Let 4> = r 1 so that 
dt = 4>'(u) du and 



4>{u) — t corresponds to r(£) = u, where r is given by 
4f(u) = l/r'mu)) 



. Now, 



From the relation, u = t(0(m)), we have that u = [<jy /4k] [e K< ^ M ) — 1] and so 



Therefore, by (jXTj) . 



V 



4k 

1 H 7T U 



V t dt 



exp(—K(f)(u))X(u)<fr(u)) du 

2 



TL 

TR A 

X(u)-^ 



TL 
TR 



V 



1 



4k 



du 



X (u)wq(u) du. 
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□ 

Although the BESQ X process and corresponding bridge are both temporally translation invariant, 
the function wq in (|5.2[) is not. Therefore, when transforming the r-space integral to a standard 

interval, "[0,r]", the integrand in (|5.ip will be modified, resulting in ("=" denoting equality in 
distribution) 



V t dt = / X(u + t l )w(u) du 
Jo 

w(u) := ai [bi + c\u]~ 2 ; a-y := 4a 



X(u)w(u) du, t = t r -t l , 



V ' 



61 := 1 + 4na v tl, c\ 



4kg, 



(5.3) 
(5.4) 



where we have suppressed the dependence on L, R in the notation, b\ and r. It is implicit, for the 
rightmost integral in f|5 .3[) . that the end values of X (at times and r) are the original values, shifted 
from times tl and tr, respectively. 

In the implementation of stopping criteria, the conditional variance of the integral in (|5.3I) is required, 
conditional on the endpoint values of the bridge. For that, we state a result on the first two moments 
of the integral of the weighted, squared Bessel bridge process. The proof of this theorem is lengthy 
and hence deferred to Appendix B. 

Theorem 5.2 For the BESQ X process X of order v, frozen on endpoints, X(0) = x and X(t) = y 
(where r > is arbitrary), and w as in i5.4]) . we have 



X{u)w{u)du\X{T) = y 



{Ax+B l )[v + 1 + -R v {- 



{B 1 +C 1 )x + {2A 1 +B 1 )y 
2r 



(5.5) 



and 

E, 



X(u)w(u)duj \X(r)=y 

= 2[A\ + B\ + A X B X -A 2 - B 2 ] + (Aj + B\ + 2(A X + B{) 2 - 2A 2 - 2B 2 )v + (At + Bi) V 
[(B l +C 1 )x + (2A 1 +B 1 )y] 2 (B 2 + C 2 - B\ )x - (3Af + B\ + 2A 1 B 1 - 2A 2 - B 2 )y 



4t 2 



(Ax + Bx) 



[(Ax + Bx) 2 + 2(A\ + B\ + A x Bx - A 2 - B 2 )]-R v 



(Bx + Cx)x + (2Ax + Bx)y , 



(Ax+Bx){ v +1 + -Rj- 



(5.6) 



where z = ^fxy, R v (r) is the Bessel quotient I u+ x(r) / I v (r), where I v is the modified Bessel function 
of the first kind (see [Tfj, ailc ^ tie constants Ai,Bi,Ci , i — 1,2, are described in Proposition \B.4l in 
Appendix B. 



Remark 5.3 It is not evident that the right-hand sides of \5. 5|) and i|5.6')) tend to zero as t tends 
to zero. However, they do; e.g., the first moment tends to zero at a linear rate and the variance 
tends to zero at a quadratic rate. To see this explicitly, we reformulate the moments in terms of 
the parameters A,b,c (see Proposition IB. 41 and Corollary IB. 31 in Appendix B) for which A —> 
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quadratically fast and c — > at a linear rate, as r — > 0. Then, in the expression for Ai,Bi,Ci 
(i = 1,2) we expand, 

2 3 

log(6 + c) = log b + log(l + c/b) = logb+ r ---^ + -^+0(T 4 ). 

Writing c = c±t, it is then straightforward to check that 

,. A x ci 
lim — = — — 

rj.0 t 2 b 2 

,. B 1 4 Cl 

l To7* = W 2 

.. Ci 2 Cl 
hm — = — — 

t|o t 2 b 2 

,. A 2 5cj 
lim — r = — 7 

r;o r 4 66 4 

S 2 8c? 



lim 

r^O T 4 156 4 

lim — r = — rr ■ 

rio r 4 36 4 



6 Adaptive Estimate of the Integral J t 3+1 V(s) ds 

We return now to the problem stated at the end of Section^ namely, the selection of the intermediate 
(quadrature) points for the estimation of the integral, J t J+1 V(s) ds. 

There are three aspects to the selection of intermediate points, at any stage in the recursion, in an 
adaptive fashion: (i) the manner of refinement (i.e., the geometric placement of an inserted partition 
point or points between those already generated); (ii) a choice of stopping criterion, to decide if 
a subinterval needs to be further refined; (iii) a decision to apply the stopping criteria locally or 
globally. By definition, a local decision means that a subinterval will be refined if the stopping 
criterion is not met on that interval, whereas a global decision means that all intervals will be refined 
if any of them do not meet the stopping criterion. 

Regarding (i), we do not force any particular refinement scheme, although for the numerical experi- 
ments in Section [51 we use bisection in t-space. Regarding (iii), we restrict our attention to the class 
of locally adaptive schemes, which we denote by ADAPT. The main purpose of the current section 
is to introduce the stopping criteria of aspect (ii). 

In general, a stopping criterion on an interval, involves a tolerance, 8, and a quantity to monitor. 
The refinement continues as long as the quantity being monitored is not within the tolerance. In 
most cases of interest, the tolerance depends on the interval being considered; in fact, it may depend 
on the entire history of refinement that led to that interval. To clarify the tolerance's dependence 
on the interval, we make a brief digression on the aspect of tolerance in a refinement scheme. 

The initial tolerance, i5o, for the "root" interval [rj,Tj + i] = [r(tj), r(tj+i)] will be user-given. For the 
sake of simplicity we shall relabel this interval as [tq, Ti] During the refining of [to, t{\, Sq is appor- 
tioned among the subintervals created by the refinement, leading to the S for each such subinterval. 
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We next describe one possible and simple set of rules of apportionment; other rules are certainly 
permitted. In particular, a more complex and efficient set of rules is described in Section [7] 

Denote the current subinterval being monitored, by [t£,tr] ee r(ifl)]. If the stopping criterion 

fails on this interval, then it is partitioned into precisely two subintervals. The typical case is where 
refinement is by bisection (either in t- or r-space) and the apportionment of tolerance is into equal 
parts. For example, if the bisection occurs in r-space, then the interval [tl, tr] inherits the tolerance, 
So(tr — 7x)/(n — r ). Similarly, if the bisection is carried out in t-space, the tolerance is given by 

So(tR-t L )/(ti-to). 

The stopping criterion for our adaptive scheme, hereafter referred to as ADAPT, is based on the 
computation of the variance of the integral. If the variance of the integral over any interval, as 
computed from (|5.3p . (|5.5[) and (|5.6|) . is below the tolerance, 5 > 0, available for the interval, we 
refrain from further partitioning of that interval: for the interval, [tl , Tr] , 



ADAPT stopping criterion : 



Var 



I V{s)du\X TL = x L ,X TR = x R ee Var / X(u)w (u) du \ X TL = x L , X TR = x,, : 

Jt L J LJt l 



< S. 
(6.1) 



This criterion is applied to each subinterval [tl,tr] and, if met, no further partitioning of that 
interval is done. If the criterion is not met, then the interval is partitioned. 

Alternatively, when the criterion is not met, the interval [tL, tR] can be partitioned; then the resulting 
subintervals can be mapped to the corresponding subintervals of [tl , tr] , for the generation of 
intermediate X values. 



7 A more efficient adaptive scheme 

In the previous section, we introduced some locally adaptive schemes. These schemes are robust in 
the sense that they recursively refine until the monitored quantity is below the relevant tolerance. 
The quantity that is monitored is subadditive in the sense that once the adaptivity algorithm has 
terminated, the sum of the contributions over all the subintervals will be smaller than the initial 
user-given tolerance, So. However, we did not address the issue that this total over all the elements 
may be significantly lower than <5o- In other words, we did not care about the minimality of our 
refinement. Due to the inherent randomness in the quantity monitored and thus the placement of 
the intermediate points, it may not be possible to obtain a minimal grid that barely satisfies the 
adaptivity criterion (I6.1[) . for example. 

However, one can do slightly better than the plain adaptivity schemes described in the previous 
section. One can approach, what we choose to describe as a near- minimal grid by introducing a global 
reservoir of tolerance, Sn- In essence, we allow the cross-subordination of the tolerance assigned to 
various intermediate elements. When an element passes the adaptivity criteria, it releases the excess 
tolerance that it has over the monitored quantity, to the reservoir. In the testing of subsequent 
elements, the monitored quantity is compared against a more lenient tolerance - one that is the sum 
of the inherited tolerance level, 5 (as usual), plus the reservoir tolerance, Sn- 

Here is the precise scheme, in algorithmic form. Suppose we need to evaluate the integral of the 
Variance process on an interval [io,ii] and further suppose the tolerance level is set to 5q. Our 
algorithm consists of the following steps: 
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Step 1: (Initialisation) Create a reserve tolerance, 5n, (whose role will be explained in Steps 3, 4) 
and initialise it to zero (5n = 0). Prepare an empty stack of 5-tuples (Tl, Tr, Xl, Xr, 6) and 
push the initial data, (t ,ti,X(T(t )),X(T(ti)),5o), onto the stack, where X is the BESQ X 
process. 

Step 2: If the stack is empty go to Step 6. 

Step 3: Pop the top element from the stack. Extract the values and check if A < 0, where 



If satisfied, go to Step 4; otherwise go to Step 5. (Notice the way in which we use the reserve 
tolerance 5-jz to facilitate passing of the stopping criterion for other "more needy" elements). 
Step 4: Set 8 n = -A. Go to Step 2. 

Step 5: Compute the midpoint T M = (T L +T R )/2. Sample X M = X(t(T m )) from the BESQ X 
process. Push the elements (Tl,Tm,Xl,Xm,S/2) and (Tm,Tr,Xm,Xr,S/2) onto the stack 
and go to Step 3. (Notice the way in which we distribute the tolerance in equal parts to the 
newly spawned intervals). 

Step 6: Exit. 

As can be seen from this algorithm, we have incorporated the following features that allow us to 
claim that the algorithm will result in a "near minimum" number of degrees of freedom necessary: 

• Our strategy is adaptive. An interval only gets subdivided if it fails the test; in the event of a 
failure, it gives out its acquired delta (from its parent interval) to the newly spawned intervals. 
Thus there is no loss (or waste) of S. 

• In the event the test is successful on any interval, it only consumes the amount of acquired 
delta that is necessary to pass the test, and releases the excess to the tolerance reservoir. The 
accumulated tolerance in the reservoir can be used later to help pass the test on the remaining 
intervals. 

Due to the inherent randomness in the above strategy, an adaptive algorithm that passes under an 
absolute minimum number of degrees of freedom, may be difficult to find. 



8 Numerical Experiments 

In this section we support the ideas introduced in this paper by a series of numerical experiments. 
To compare the results of our method with other methods available, such as the finite difference 
method, we choose the same test problem as the one appearing in [5] (Table I, Section 4). For ease 
of reference, the input parameter values are: S = 100, K = 100, V = 0.010201, re = 6.21, = 0.019, 
ay = 0.61, p = —0.7, r = 3.19%, T = 1.0 year. Also, as a benchmark, the true option price = 6.8061. 

In our tests we compare our methodology with a variant of the finite difference method which uses 
a predictor-corrector step for better convergence. Also, we have chosen to carry out the refinement 
scheme in the t-space and the refinement involves simple bisection of the interval in question. 

Since our method is based on discretisation, a bias in the numerical method is expected. However, 
in contrast to the other methods that we know of, our method allows a systematic control of this 
bias by means of relating the bias to another numerically observable quantity: the variance of the 
integral of the variance process over the interval in question. Furthermore, we put great effort in 
ensuring that the computational cost of achieving the limit on the variance of the integral (and 



A := Var 




V( S )ds\X(T(T L )) =X l ,X(t(T r )) =X r 



(5 + 5n) 
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therefore indirectly on the bias in our method) is kept to a near minimum as explained in Section [7] 
To this end, it is clear that the very first set of tests should demonstrate the ability of our algorithm 
to consume as little computational resources as possible. We do this in Tables Q] and [2] For our 
test problem, we note that the dimension of the BESQ process is given by A = 1.2684 and the 
order is v — —0.3658. The origin is accessible and is also reflective. The moderately high value of 
the volatility of the Variance may also result in many paths taking very high values. Based on this 
observation we decided to capture the properties of our algorithm by banding the endpoint value. As 
expected, we see from Table [1] that for a path that starts at a moderate value and reaches a fairly 
high value, the number of the intermediate points that need to be inserted, to reduce the variance 
of the integral below a given tolerance, is also very high. It appears that the number of intermediate 
points is proportional to the absolute difference between the left and the right endpoint values. We 
also observe that the distribution of the number of intermediate points is tighter for lower endpoint 
values compared to higher endpoint values. 

In Table [5J we demonstrate the effect of cross subordination of tolerance which forms an essential 
part of our algorithm. What we show in Table [H is the tolerance level that was wasted by the 
algorithm; i.e., for the interval in question, the difference between the given tolerance by the user, 
and the sum of the actual variances of the integrals over the subintervals. We again study this in 
terms of bands of endpoint values. As to be expected, the wasted tolerance decreases significantly 
when more and more intermediate points are inserted (as can be seen with the distribution's very 
short left tail, for high endpoint values). This effect is clear because, when inserting more points, 
more iterations are made and hence more use is made of the reserve tolerance. This is an important 
feature of our algorithm as it shows that the wastage is minimal when there is highest demand for 
refinement. 

Figures [1] and [2] illustrate the expected fact that the bias is reduced in both the predictor-corrector 
method and our adaptive algorithm, as more and more intermediate points are introduced. In the 
predictor-corrector method the independent variable is directly the number of intermediate parti- 
tions, whereas in our algorithm, the independent variable is the tolerance provided by the user. Due 
to this mismatch in the independent variable it is not obvious how to compare the relative perfor- 
mance of these methods. So in this sense, Figures [T] and [2] may be considered just a sanity check of 
the expected way in which these algorithms are supposed to work. 

Our next task is to compare the predictor-corrector method with our algorithm. For this we in- 
troduce the quantification of accuracy, which is defined as the inverse of the absolute relative bias. 
Furthermore, we saw from Figures [1] and [5] that using the predictor-corrector method with smaller 
interval size has the same directional effect as reducing the tolerance level in our algorithm, and 
both these actions result in increasing the time spent in simulation. Therefore the most ideal way 
of comparing the two methods is by plotting the accuracy, as defined above, versus the time spent 
in simulation. This is what is shown in Figure [31 As a byproduct of this analysis, we make a very 
interesting observation regarding our method. It is exponentially rewarding in the initial part with 
a much steeper slope than the predictor-corrector method. Note that the accuracy is plotted on 
a logarithmic scale. We also note that both the methods taper off as we move to the right, with 
diminishing rewards. This may be attributed to the fact that reducing the bias substantially below 
the simulation error is fruitless. This last observation brings us to the final figure of this section, 
Figure HJ In most risk management work, the computational budget associated with pricing a fi- 
nancial derivative is limited. Based on that constraint, only a small number of simulation paths 
(typically between 1000-5000) are used for pricing. Of course, this results in a large value for the 
error associated with the Monte Carlo method. It is clearly pointless to control the bias to any order 
of magnitude below this error. Figure [4] allows us to demonstrate that controlling the bias with a 
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Table 1: Number of intervals used in the calculation of the integral of V (tolerance = 0.000001). 



right bin 


Endpoint Variance 


boundaries 


0.000001 


0.0001 
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0.09 
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385 


594 


1016 


1658 
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3037 


3027 
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1057 


56 


12 


19 


36 
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244 


463 


908 
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2279 


2801 


2999 


2453 


64 


3 
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22 


70 


141 


298 


646 


1124 


1744 


2342 


2820 


72 





1 





2 


8 


33 


85 


203 


419 


726 


1266 


1973 


80 











2 


1 


6 


21 


59 


139 


293 


536 


1016 


88 














2 


3 


10 


41 


65 


128 


223 


518 



Table 2: Unused tolerance (0.000001) in the calculation of the integral of V. 
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boundaries 


0.000001 
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80 
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59 
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2 


3 


10 
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65 
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223 
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tolerance below 1.56c-06, say, has diminishing returns. 
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Figure 1: Estimation of bias in the predictor-corrector method (1000 trials of 10000 samples). 
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Figure 2: Estimation of bias in the ADAPT method (1000 trials of 10000 samples). 
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y- Predictor Corrector ^-Adaptive Bias Control 



~ay 100 




Figure 3: Comparison of the predictor-corrector against ADAPT algorithms: accuracy versus time; 
logarithmic scale on vertical axis. 



(Based on 10000 Trials of 1000 Samples) 




Figure 4: Accuracy versus time trade-off in the ADAPT method (based on 10000 trials with 1000 
samples). Node labels on graph are the tolerances logarithmic scale on vertical axis. 
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9 Conclusion 

In this paper we explored a variant of the methodology proposed by Broadie and Kaya in [5], to 
simulate the dynamics of the Heston model. As our method relies on the numerical computation of 
the integral of the variance process, it is subject to bias. The adaptive nature of our method allows 
the efficient, practical control of this bias. 

It is expected that for options on instruments with a large number of reset dates, the adaptive 
method of this paper should outperform the exact method of Broadie and Kaya. 

In summary, our method provides the following benefits: 

1. Unlike finite-difference methods, our method cannot generate negative values for the Bessel 
process and associated bridge. 

2. Bias is efficiently controlled, as shown numerically by Figures [5] and @] 

3. It allows a much greater degree of flexibility than any of the other methods we have seen. This 
is advantageous where one can increase the tolerance level for middle-office risk-management 
work and reduce the tolerance level for front-office pricing. In other words, the tolerance level 
is a function of computational budget that one has at one's disposal. 

4. Perhaps the most interesting feature of our method is a "near-invariance" of the number of par- 
titions (intermediate reset /coupon dates). Based on the observation that the most demanding 
part of our algorithm is the adaptive computation of the integral, we expect that the number 
of intermediate points we require for one big step of T years is roughly equal to the total 
number of intermediate points required, had we decided to take n steps of length T/n. In the 
latter case our method should perform much better than the Broadie-Kaya method. 

5. As our algorithm derives, in essence, from the Broadie-Kaya algorithm we can safely assume 
that all the extensions to jump diffusion models presented in [5] should work with our algorithm 
as well. For the sake of brevity we do not reproduce the details in our paper. 

6. Our method is straightforward to implement because many of the generators are now readily 
available in standard libraries. 

The most important extension that awaits investigation, is to higher dimensional systems, for pricing 
options on equity baskets. 



Appendices 

Appendix A: Proof of Corollary 13.21 

We will actually give the derivation in the opposite direction, from X to V. Let <f> = t _1 so that 
u = r(t) corresponds to t = 4>{u)\ uq = r(to). Also, for a constant c > 0, to be determined, denote 
Y t c = e ct X T(t) = e^Xu, so that Y t ° = X u . Now, by Q , 




so 
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where we have used a general time-substitution result for stochastic integrals (see Proposition (30.10) 
in Chapter IV of [2]) to transform the integral in the second equality. Therefore 



dY t c 



cY t c dt + e ct [Xdr{t) + 2^X r{t) dW r{t) ] 



cY t c dt + X^fdt + 2e ct/2 JYfdW T 



(*)■ 



Let Z be a given Brownian motion and take for W, the process defined by 



2 



so that 



and, by Ito's formula, 



J<j>(u ) 

W r{t) = ^ fe-^dZ s 

1 Jta 



dW 2 T{t) = 



W T{t) e- ct/2 dZ t 



dt. 



Thus 



W 2 v 



dt = W 2 (t) - -f[exp(-ct) - exp(-ct )]/c 



(A.l) 



(A.2) 



is a martingale. With the choice c = — k, we obtain that W 2 ^ — r(t) is a martingale and changing 
variables back to u, that W 2 — u is a martingale (with respect to a different filtration, of course). 
Since W itself is clearly a continuous martingale, we conclude from Levy's theorem that W is a 
Brownian motion on [to, oo). 

Returning to (|A.1|) and substituting the differential form of (IA.2I) . we obtain, with V = Y~ K and 



dV t 



-KV t dt + n6dt + <J V \/Vt dZ t 



Finally, we identify Z with 



□ 



Appendix B: Proof of Theorem 15.21 

Our starting point is the following representation taken from (T3] (see Theorem 3.2 and its proof 
on pages 442-443, therein^) for the Laplace functional of the integral of a BESQ bridge, X, which 
starts at x and ends at y: 



E, 



exp 



X(u) dfi(u) 



E, 



exp 



X(u)dp(u)} X(l) = y 



x,y>0 



where /x is a Radon measure on [0, oo) with support in [0,1]. 



2. Note that we have replaced the arbitrary measure /i by 2/i and corrected two typographic errors on page 443: 8 
should be divided by 2 and the subscript p 2 (l) on q, should be tr 2 (l). 
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Theorem B.l Let /.i be a Radon measure on [0, oo) with support in [0, 1], and set 



F(x,y,n) := E xy 



exp 



X(u) d/j,(u) 



Then 



F(x,y,n) 



■ exp< 



1) / (/){u)- 2 du 



1 /„ 



0'(O) - (J 4>{u)- 2 du\ +1 



• exp< 



l - -I' 

-2 



1 - <j>(iy / <j)(u)- z du 



where I u is the modified Bessel function of order v and is the unique solution (in the sense of 
generalised functions) of the ODE 



2/i ■ 0, (f)(0) = 1, > 0, nonincreasing on [0, oo). 



Consequently, is convex and right-differentiable with a right-continuous, nonpositive right-derivative. 



This result can be transferred to general interval [0,r], by using the following result (see [T2] or 
[T5]): If QJJ y denotes the law of X under which it is a Bessel bridge, starting at x and ending at y 
(at time r), then y is also the <3^/ T y / T -law of tX(u/t), on [0, r]. The only case of interest to us, 
is when // has a density with respect to Lebesgue measure, necessarily of the form, m(u) lro iT i(tt), 
< m < oo. With a slight abuse of notation, we write F(x, y, m) instead of F{x, y, fx) in this setting. 



The details of the transformation are as follows. Given the density m, we set m T (u) 



m(ru) 



< u < 1; m T (u) := 0, for u > 1. Using the cited equivalence of laws and then making the change 
of variables, u i-> tu, yields (with W xy denoting the Ql y expectation): 



Ky 



exp ^ — j X(u)m(u)du 



= E 



x/t,ij/t 



E l/r,y/r 



exp 



F (-, -,m T 

\T T 



exp < — / TX(u/T)m(u) du 



X(u)m T (u) du 



There is one more technical result which is needed to completely localise the problem to the support 
of /i. Again we restrict attention to the case where fi has a density m which we assume is continuous 
on its support, the interval [0, 1]. In that case, the ODE is 

(j)" = 2mlfo j i] • 4>, 0(0) = 1, 4> > 0, 4> nonincreasing. 

Standard regularity theory yields that 4> is smooth (C 2 ) on [0, 1) and satisfies the equation, <\>" = 2m0 
in the classical sense thereon. Of course, <fi is constant on (l,oo) and being continuous everywhere, 
the constant value is 0(1). We now show that satisfies a Neumann boundary condition at u = 1. 

Lemma B.2 The left-hand derivative, 0'(1~) = 0, so that is C 1 smooth across u = 1. 



1G 



Adaptive Simulation 



Proof Let g 6 Cg°((0, oo)) be a test function. The ODE (aside from the boundary and side 
conditions) means that 



Also, 



4>{u)g (u) du = / m(u)l[o,i] (u)(f>(u)g(u) du = I m{u)(j>{u)g[u) du 



1 />oo />1 />oo 

- / 4>{u)g"{u)du = / ^(u)g"(u)du+ / 0(w) 5 » 
* Jo Jo Jl 



(u)g"(u)du + (f>{l) / 



<P(u)g"(u)du-ct>(l)g'(l) 



= 4>(l)g'(l)- ^(u)g'(u)du-cf>(l)g'(l) 



= -^(l-)fl(l) + / <l>"(u)g(u)du 



</>'(! ).g(l) + / m(u)0(u) fl (u)du. 



Therefore 0'(l~)g(l) = for all g <E C§°((0, oo)), which implies that <f>'(l~) = 0. 
Thus we arrive at the final formulation of the required Laplace transform: 



□ 



Corollary B.3 Let m T = a[b + cu] 2 , < u < 1, where 



88t 2 



, b=l 



(m T = 29t 2 w(tu), where w was introduced at (15.41) ) and set 



Then 



L(0) := E x 



exp < — 8 / X(u)w(u)du 



L(9) 



I" 1 h 



(/>(!) / cj)(uy z du 



-0(i)J o V(n)-^~ 



• exp< 



• exp< 



x 

JL 

2r 



0'(O)- / 4>{uy 2 du) +1 



i - 0(i) 2 / (f>(uy 2 du 



(B.l) 
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where I v is the modified Bessel function of order v and <fi is the unique solution to the BVP: 

<j)" = m T {u) ■ 0, 0(0) = 1, 0'(1) = 0. (B.2) 



For the first two moments, we are interested in the coefficients of —6 and 8 2 /2 in the Taylor expansion 
of L(8) about 8 — 0, since e~ r = 1 — r + r 2 /2 + • • • and thus the left-hand side of (|B.1[) equals 



1 - 6 ■ E„ 



X(u)w(u) du 



X(u)w(u) du 



0{8 3 



Accordingly, we work out the Taylor expansion of the right-hand side of (jB.ll) to second order. 

With reference to (|B.2j) . the function depends implicitly on 8 through the constant, a, which 
appears in the definition of the function m T . We make this dependence explicit in our notation, by 
writing cj)(u;8). Clearly 0(- ;0) = 1; so <//(0;0) = 0. (Differentiation with respect to u will continue 
to be denoted by a prime (') superscript; differentiation with respect to 8 will be written explicitly; 
e.g., as a partial derivative, d/d8.) Thus we set 

0(1; 8) = l + A 1 8 + A 2 8 2 +0(8 3 ) 

(u;8)- 2 du = 1 + Bx6 + B 2 8 2 + 0{8 3 ) 

0'(O;0) = d8 + C 2 8 2 + 0{8 3 ) 

leaving the determination of the coefficients, Aj, Bi, Ci (i=l,2) for later. 

Expansion of the right-hand side of ([B.ljl in powers of 8, can be effected in a few stages. The terms 
involving a multiplicative inverse, like (J Q 4>{u)~ 2 du)~ l , can be expanded using the expansion for 
and the geometric series expansion, = (1 — [1 — r]) _1 = 1 + [1 — r] + [1 — r] 2 + 0([1 — r] 3 ), for 
r close to 1. The two exponentials can be combined and then handled with the usual expansion, 
e r = 1 + r + r 2 + 0(r 3 ), for r close to 0. The ratio of modified Bessel functions can be handled, using 
the following two identities for modified Bessel functions (see 9.6.1, 9.6.26 in pQ): 

,2\ 1 

(B.3) 
(B.4) 



I'M = 1 + -\l v {r)--I' v {r) 
\ r A J r 

I' v {r) = I v+1 (r) + ~I«(r). 



We can substitute (|B .4|1 into (|B.3I) and divide both identities by I v {r) to obtain the following ones 
in terms of the so-called Bessel quotient function, R u = /„ + i(r)//„(r): 



= 1 



--Mr) 



- + R v (r). 



We can apply these identities to (|B.1[) by writing r = ifxy/r and r = r o /0(l) J Q 0(u) 2 du, and 
expressing the ratio of Bessel functions in (|B.1|) , in the form 

I»(r) _ J„(r + [r - r ]) _ I v {r Q ) + I'M[r - r ]) + I'J(r )[r ~ r Q } 2 /2 + 0([r - r ] 3 ) 



I v {r Q ) 



I v {r ) 
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and noting that 

r-r = r [(6(l) f 6{u)- 2 du)- 1 - 1] 
Jo 

has an expansion in 0, without constant term. 

The remaining details are straightforward but tedious algebra which is omitted. The final result is 
the expression in Theorem 15.21 

We now turn to the calculation of the constants, Ai, Bi, Ci (i=l,2), in terms of the constants a, b, c. 
Denoting 

86(u;0) 



Ag(u) := 

r*(u) := 



80 ' 
8 2 6(u;0) 
80 2 



A{u) := A (u) 
, T(u) := T (u), 



we then have 



A, = A(l), A 2 = ir(l); d = A'(0), C 2 = lr'(o). 



Also, for Bi,B 2 , note that 
— t 

dej 

d0 2 J o 



u;0) 2 du = 
u\0)~ 2 du = 



&<j>(u;0)- 



86 
86{u; 0) 
80 



(B.5) 



Evaluating these results at = + , we obtain 



Bl d0 



(u-0)~ 2 du 



B 2 = 



6=0+ JO 

2d0 2 



6=0+ ^0 



= -2 / A(u)du, 
Jo 

du = [ 3A(u) 2 - T(u)du. 
Jo 



(B.6) 
(B.7) 



The functions A and T can be found by differentiating {-§q) the ODE and boundary conditions for 
6 and then solving the resulting, very simple problems at = + . To that end, we bring out the 
^-dependence of m T explicitly by writing 

m T = 0M T 

with M T independent of the parameter, 0. Then, 

8 8 8d) 

K = = qq{0M t (I>) = M T 6 + 0M T -^ =^ A" = M T ; also A(0) - 0, A'(l) = 0; 



8 2 8 6 8 2 S" 



T" = 2M T A; also T(0) = 0, r'(l) = 0. 



Integrating and using the boundary conditions, we obtain: 

A'(u) = - [ M T (v)dv, A(u) = I A'(v)dv 



T'(u) = -2 / M T (v)A(v)dv, T(u) = / T'(v)dv. 



(B.8) 
(B.9) 
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Proposition B.4 With a,b,c as in Corollarv \B.3[ A := a/8, and log 2 denoting the square of the 
log [unction, 



Ax 
A 2 
Bi 
B 2 

Ci 

c 2 



A 

c 



' +-\og[b/(b + c)] 



b + c c 
A 2 



c 4 (6 + c) 
2A 



{b + c) log" [6/ (6 + c)] - 361og[6/(6 + c)] 



+ log[6/(6 + c)] 



c(3& + 2c) 



2(6 + c) c 



(B.10) 

(B.H) 
(B.12) 



-A 2 



c 5 (6 + c) 

• [2(6 + c) 2 log 2 [6/(6 + c)] - 2(36 2 + 8c6 + 4c 2 ) log[6/(6 + c)] - c(66 + 11c)] 



6(6 + c) 
A 2 
c 3 (6 + c) 



c(26 + c) 
6(6 + c) 



+ 21og[6/(6 + c)] 



(B.13) 
(B.14) 

(B.15) 



Proof All of the integrals in (|B.8I) and (|B.9I) are straightforward to evaluate, as well as the integrals 
(|B.6|) and (|B.7p , for i?i and -82- We omit the elementary calculus and just state the end results, as 
they can be easily verified by differentiation and checking a boundary condition. Then one can use 
(|B.5P to obtain Ai,A 2 ,Ci and C 2 . 



A'(«) 
A(u) 

t 

A(v) dv 

r» 

T(u) 



- [(6 + c)- 1 -(6 + cu)- 1 ] 

a r 1 



(6 + c) 1 u+ -log[6/(6 + cw)] 



2(6 + c 
2A 2 



+ 



(1 + log b)u 6 + cu 



log [6 + cu] + — log 6 



c 3 (6 + c) 



2A 2 
c 3 (6 + c) 





' 6 " 


log 






6 + c_ 



6 + c 
6 + cu 
6 + c 



log 



log[6 + c] 
6 



c^ c^ 
26 + c 26 + c 



2c 
6 + c 
~ 2c 
6 + cu 



(log 6)^ 



6 + cu 
36 + c 



log 6 



6 + c 6 + cu 
log [6 + cu 

26 + c - (6 + c) log 6 



log [6 



(log [6 + cu]) 2 + log 



36 + 2c 



'(6 + c) 2 



log [6 + cu] 
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3A(w) 2 - T(v) dv 



A 



2„,3 



A 2 



c 2 (b + c) 2 cJ>(b + cf 



[2c 3 + 36c 2 - (6 + c)c 2 log b + 2(6 + c)c 2 log [6 + c]] i 



2A 2 6(36 + 2c) log 6 ^ 2 [-4(6 + c)6 2 + (56 + 3c)6 2 log 6 - 66c(6 + c)] 



c 4 (6 + c) 2 
A 2 



c 5 (6 + c) 

A 2 
c 5 (6 + c) 



c 5 (6 + c) 2 

(6 + cu) 2 [3 log[6/6 + cu] - log[6 + cu] + 2] 



2 (6 + cu) [2(6 + c) 2 log 2 [6/6 + cu] - 2 (6 2 - 3c6 - 3c 2 ) log[6/6 + cu] 
- 26(26 + c) log [6 + cm] + 2(6 + 3c) (6 + c)] . 



□ 



Appendix C: Explicit Laplace transform 

In the previous appendix, we avoided solving the BVP (|B.2|) for </> and then calculating J Q 0(tt) -2 du. 
In order to extricate the moments of fT X(u)w(u) du from its Laplace transform, it was sufficient 
to apply the method of variation of parameters, which led to expressions for the moments in terms 
of 4> and its derivatives with respect to the Laplace parameter, 9, at 9 = 0. This method could work 
for more general weight functions than our specific w. 

However, to distributionally validate the ADAPT schemes, it is most convenient to have an explicit 
form for the Laplace transform (|B.1[) . and for that we now solve the BVP. 

Proposition C.l The solution to the BVP 

</>" = m T (u) ■ <j>, 0(0) = 1, <t>'{l) = (C.l) 



is 



(j)(u) = eb p (b + cu)- p + (l-e)b q (b + cu)- q (C.2) 



-1 + \Jc 2 +~4a/c -1 - Vc 2 +4a/c n , 
P = 5 . <1 = o ( C - 3 ) 



e = -qb q /[pb p {b + c)-^^ -qb q }. (C.4) 



where a, 6, c were defined in Corollarv \B.3l 

Proof It is a simple matter to verify that the function described by (|C.2[) - (|C.4[) satisfies the 
ODE and boundary conditions (jC.ip . A sketch of the method of solution is as follows. We set 
ip(u) = <j)(u/yfa) and then tp(u) = f(r) where r = [6 + <5tt] _1 and 6 = c/yfa. This leads to the 
following BVP for / on the interval [r e , 6 _1 ], r e = [6 + c] _1 : 
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Seeking a solution of the ODE alone, in the form r p , leads to the condition 

p 2 +P-5-' 2 = 

which has the two solutions (|C.3I) . The function 

/ (r) = eb p r p + (1 - e)b q r q 

then satisfies the boundary condition, = 1, for any e, while the choice, (|C.4|) . guarantees that 

the other boundary condition, f'(r e ) — 0, is satisfied. Unwinding the definitions from / back to cj> 



yields (fCT2|> . □ 
Corollary C.2 With p, q defined at if03)) . 



Jo 



<j>'(0) = -f- ( (C5) 

be p bP(b + c)-^+^/ c - qbi 

m = I P(b + c)- p -qib + er q 

b pbP(b + c)-^+^/ c -qbi 



-2 



du 



(b(b + c))i-Pc 2 



8 (c 2 + 4a6>) 3/2 ((c 2 + 2a<9) (6(6 + c))*-* + a (6 2 (p-?) + (6 + c) 2 (J>-9)) 0) 

• f(p - g)6 p -«+^ - 6f-?+^ + (b + c) p - q Vb +(b + c) p - q Vb{p - 9 )) 2 

• (2c 2 (p - g)(6(6 + c)) p -« + (b + c) 2 ( p - q ) (c 2 + 4ad) 

- b 2(p - q) {Aa6 + c (c + c(p - ?))) - c(6 + c) 2 ( p -^c(p - q)) (C.7) 



Proof The results (|C5|> and (|C6[) are straightforward calculations. On the other hand, the deriva- 
tion of (|C.7j) is not straightforward; it was effected with the aid of Mathematica® . □ 
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